Operator folding and matrix product states in linearly-coupled bosonic arrays 
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A protocol to obtain the matrix product state representation of a class of boson states is in- 
troduced. The proposal is presented in the context of linear systems and is tested by performing 
simulations of a reference model. The method can be applied regardless of the details of the coupling 
among modes and can be used to extract the most significant contribution of the tensorial repre- 
sentation. Characteristic issues as well as potential variants of the proposed protocol are discussed. 



PACS numbers: 

INTRODUCTION 

The realization that quantum states can be written in 
terms of a tensor network whose elements display inter- 
esting properties has prompted a wealth of research in 
what is nowadays known as the field of Matrix Product 
State (MPS) [if. Although the properties of MPS can 
be exploited in a variety of ways, it is Time Evolving 
Block Decimation (TEBD) and Density Matrix Renor- 
malization Group (DMRG), together with their variants, 
that have proved highly robust and appropriate in most 
situations of interest. However, other methods have also 
been proposed, for example, in the area of infinite chains, 
where the calculation of local mean values can be formu- 
lated in terms of bundled tensor networks, or in the area 
of Gaussian states, where the MPS network is obtained 
as projections of highly entangled states [2|. MPS of- 
fers a view that is particularly convenient in variational 
approaches, where some physical state is obtained by 
renormalizing a tensor network. This has led to an inter- 
est in classes of states that can be efficiently simulated 
[3|. Notwithstanding its recurrent use in spin models, 
the relevance of MPS is especially notorious in bosonic 
systems. In this context, the application of TEBD has 
allowed the numerical exploration of boson chains under 
different conditions [^-[Sl revealing phases and regimes 
with very interesting properties. 

Perhaps the most elementary way of representing a 
quantum state is as a set of complex coefficients derived 
by writing such a state as a superposition of elements 
of a basis. In what respects to indistinguishable parti- 
cles, the basis is constituted by occupation states upon 
which ladder operators can raise or lower the associated 
number of particles. Because any of these states can be 
put in terms of ladder operators acting on the vacuum, 
it is possible to envisage a representation relative to such 
operators. This approach is practical, for example, when 
the symmetries of the problem allow an advantageous 
handling of the Heisenberg equations [9[ . This is seen in 
linear systems where the underlying physics is driven by 
interference and single body (SB) effects. These systems 
are quite recurrent, not only as realistic descriptions of 



physical phenomena, such as optical fields [10[ or weakly- 
interacting Bose-Einstein condensates, but also as mod- 
eling tools. The latter case is manifest, for instance, in 
the framework of the mean field or Hartree-Fock approx- 
imation. Insight in this direction must therefore be of 
significance 

In the development that follows, a method is proposed 
to go from a representation of a bosonic state in terms of 
operators to a canonical MPS representation. The anal- 
ysis makes use of the properties of both representations 
and the central argument does not involve approxima- 
tions. Results obtained using the proposed technique are 
compared against benchmark data. It is pointed out that 
the range of applicability does not depend on boundary 
conditions or number of next-neighbors, but rather on 
whether the state can be put in a compatible form. In 
the final part, potential applications and complementary 
remarks are set forth. 



LINEAR BOSONIC SYSTEMS 

Following a second quantization scheme, let us propose 
a system of M bosons. Every boson can occupy N quan- 
tum levels which are characterized by the bosonic opera- 
tors hj and a\ satisfying [aj,a{.] = 6^; and [aj,ak\ = 
with j,k = 1,2,...,A^. In absence of interaction, the 
Hamiltonian can be written as 



N N 



H ^^^ hj.ka]dk, hj,k = /ifcj. 
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Matrix hj^k [h) is the Hamiltonian when I\I = \. h 
also defines the operator dynamics according to 
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which can be obtained by differentiation of a- = 
;^'*^a^e**-^ {h = 1). A product of local Fock states 
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where |0) is the state with no bosons. More com- 
plex configurations can be constructed as superposition 
of these states. Now let e/ be an eigenvalue of h corre- 
sponding to the normalized eigenstate |e/) 
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hj,k£k,i = eiejA {1,3 = l,2,...,iV). 



(4) 



An eigenstate \Eni. 
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of H with cigencnergy 
+ UNf-N can be built as a 



product of SB cigcnmodes as 
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(5) 



The size of the basis is {N + M - iy./M\{N - 1)!. It 
can be seen that the state of a system of free bosons is 
determined fundamentally by the contribution of the SB 
Hamiltonian and the interference effects arising from in- 
distinguishability, which is implicit in the bosonic opera- 
tors. This characteristic renders the system into a linear 
regime, where a composition of solutions of H, like in 
Eq. dl]), is also a solution, and a SB eigenmode remains 
physically unaffected by other SB eigenmodes. 



BOSONIC STATES IN MPS FORM 

In order to establish a ground to perform the transition 
to MPS, let us imagine that bosons are arranged in a 
chain with open boundary conditions. This assumption 
however does not need to coincide with the real boundary 
conditions of the problem. A site in the chain is labeled 
by the integer n ranging from 1 in the right end to N 
in the left end. Using MPS, the quantum state can be 
represented as a superposition of non-local states in the 
following way (up to few changes, the notation in ll| is 
followed) 



|^) = ^AL"lrNb)A{r 



-l]|^[A^:"+l])|p["])|^["-l:l] 



(6) 



of the chain. AJT" and AIJ* are the Schmidt coefficients 
associated to such decompositions. The states |p'"l) are 
elements of a local basis at site n. For bosons, it is con- 
venient to choose a local Fock basis. The complex coeffi- 
cient rL"/j (p) determines the contribution of a basis state 
to the superposition. Integer p is an occupation number 
and ranges from to M. Integers /i and i' are labels of 
two distinct sets of Schmidt vectors. The maximum num- 
ber of these vectors over all possible bipartite decomposi- 
tions of the chain is called x- An important aspect of the 
MPS representation is that by adjusting x it is possible 
to control the number of coefficients employed to describe 
the state. This allows to approximate huge states by re- 
taining the most significant contribution of their respec- 
tive MPS representations (the part linked to the biggest 
As). The set of tensors {rL"|i(p), A|I' .V(//, z^,p, n)} is a 
representation of {ip) that can be updated when an uni- 
tary transformation is applied on a pair of consecutive 
sites. In what follows, it is shown how this feature can 
be applied to put states like ([3]) or ([S]) in MPS form. 

Let us start by considering the simplified case where 
ni bosons occupy the same arbitrary SB state. The state 
can then be written in terms of a non-diagonal mode 
(NDM) as 



|-0) = ^== ci,ia{ + C2,ia^ + ... +CAr,ia^ |0). (7) 

The meaning of the second subscript in the coefficients 
is explained further down. Normalization of j^/') requires 
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In a first step all these coefficients are to be made real. 
The idea is to operate on | V') with a series of local unitary 
transformations that act on the operators and take away 
the complex phases of the coefficients as follows 



— i0i_ia, Oi ~tgJ0i,iaj ai 



-i4>l,l?,^ 



^Q4^|c/,i|, (9) 



where (pi^i is the phase of c;_i. This is done for / = 
l,2,...,iV. The order in which the transformations are 
applied is not important. Next, a rotation operation is 
applied on a couple of neighbor sites using the angular 
momentum operator 



|^["-i^i]) and |jy[^^"+i]) are, in that order, Schmidt 
vectors to the right and left of site n (superscripts in- 
dicate the vector subspace). Notice that on each case 
such Schmidt vectors belong to different decompositions 
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Explicitly, this transformation reads, 
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Consequently, the contribution of al+i can always be 
suppressed by choosing the appropriate angle, namely. 



tan 



2 






(12) 



If the procedure is first utilized to suppress aL, then 



^7V' 



one can successively suppress the other ladder operators 
in decreasing order until just (aj^)"! is left acting on |0). 
Here, this process is referred to as Folding. The inverse 
process, or Unfolding, is just a way of getting the original 
state back 




g J-^ + 1. 
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-|0). (13) 



Notice that now the order in which two-site transfor- 
mations are applied matters. The order of multiplication 
is assumed to be 
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which requires N' <N. To fold ([15]), the first NDM is 
folded as shown for ([7]) . This affect the coefficients of the 
other NDMs but because the transformations are linear 
in the operators, the new coefficients obey a relation like 
([TS]) . As a result, after folding the first NDM, the coef- 
ficients of a\ automatically vanish in the other NDMs. 
The process can then be applied again to fold the second 
NDM, but this time it is more reasonable to fold until 
flj is left alone, skipping the last folding operation, since 
a[ is not present in the second NDM. In this way, fold- 
ing the second NDM does not unfold the first mode and 
0.2 disappears from the rest of NDMs. The procedure is 
repeated in a similar way until the state is reduced to a 
simple product of local Fock states. The original state 
([T5|) can therefore be recovered as 



1 N N-l 
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which in turn can be numerically implemented in terms 
of MPS as explained before. 



APPLICATIONS 
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(14) 



and analogously in subsequent expressions. The sec- 
ond subscript in the angles and the coefficients makes 
reference to the only mode left after Folding. In order to 
write (|13p as a set of tensors, the state with ni bosons in 
the first place of the chain is written as MPS. This can be 
readily done because the Schmidt vectors of such a state 
have a simple structure. Subsequently, the tensorial rep- 
resentation is updated according to Eq. ([T3|) . following 
the protocols available for one- and two site operations 

m- 

More complex situations take place when bosons are 
distributed over several SB states. As has been seen, an 
important class of these states can be generically repre- 
sented as 
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together with 



In order to test Unfolding in a controlled manner, a 
Hamiltonian with a known analytical profile is brought 
up, namely 



ijM 



"k+l 



5i^\ 



(18) 



plus periodic boundary conditions, hj_N+i ~ hj^i and 

hN+i,k = hi.k (ajv+i ^ '^i)- ^^ Ih*^ spectrum of this 
Hamiltonian is in general degenerate, the next reference 
eigensystem is chosen 



e; = 2 cos 
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(19) 



The calculation consists in solving Eq. ([2]) and then 
inserting the dynamical operators in Eq. (j3|), assum- 
ing that at i = there is one boson at each site of the 
chain. The resulting state is then written as a tensor 
network using Unfolding. To do this, Eq. pT]) is im- 
plemented as a numerical routine that integrates the up- 
dating subroutine of the programs described in [8[ . The 
obtained results are then compared against equivalent 
simulations carried by diagonalization. The lower panel 




FIG. 1: Entropy between one site and the rest of the system 
(Top) and error A = 1 — \{ij}\tp')\^ (Bottom) in a boson chain 
with N = 8 and M = 8 initialized with one particle at each 
site. The underlying Hamiltonian displays next-neighbor hop- 
ping (Eqs. IT]) and (|18|) ') and the boundary conditions are pe- 
riodic. The error determines the difTerence between the state 
found by standard diagonalization (|i/'')) S'nd by Unfolding as 
explained in the text. 



of Fig. [T]shows, as a function of time, the numerical error 
produced by Unfolding when compared to the standard 
method. Unless otherwise stated, it must be assumed 
that in the MPS computations x is not bounded but dy- 
namically determined by the updating routine as the sim- 
ulation runs. In this way, all the elements of the MPS 
representation are retained. As can be seen in Fig. [1] er- 
ror is comparable to computer precision and it does not 
grow over long intervals. This because in Unfolding the 
state for a given time only depends on the initial condi- 
tion and the solution of the equations of motion for the 
operators, which can be obtained with high accuracy for 
any t. The upper panel of Fig. [1] shows the single site 
entropy of the chain, calculated from 



s--j:M i«g(Ai^' 



(20) 



S measures the entanglement between one site and the 
rest of the chain and can be easily computed from a MPS 
representation. It is known that the chain relaxes to a 
Gaussian state with maximum entropy subject to fixed 
second moments [1^. As a result, the saturation of S 
determines a time window along which the dynamics is 
relevant. 

Fig. [2] shows S for the eigenstates of H as well as the 
numerical error incurred by passing such eigenstates to 
MPS using Unfolding. In this figure every state has been 
represented only by the exponents that appear in Eq. ([5]) . 
This can be done because a SB eigenmodc formed from 
pop can be transformed into any other SB eigenmodc 
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FIG. 2: Entropy between one site and the rest of the system 
(Top) and error AE — \E — -Eni,...,7ijv I (Bottom) of energy 
eigenstates of Hamiltonian (1) for N — 8 and M — 8. The 
eigenstates are built as products of SB eigenmodes, as de- 
scribed by Eq. ([5]), and then converted to MPS in order to 
find S and E. Since S depends only on the exponents of the 
product, a many-particle state is represented by the number 
of bosons at each SB state, making no reference to which SB 
state the exponent actually apply. For instance, 17 means two 
SB states are involved, the first with 1 boson and the second 
with 7 bosons. Entropy is independent on the specific choice 
of such SB states. 



of the same family using only single-site unitary opera- 
tions. Recall that invariance under local unitary transfor- 
mations is a property of entanglement. Fig. [2] suggests 
that eigenstates of H made of bosons distributed over 
many SB eigenstates contain more entanglement than 
eigenstates with bosons arranged over few SB eigenstates. 
Nevertheless, the eigcnstate of H with all SB states oc- 
cupied docs not show maximum entanglement. 

The efficiency of Unfolding as a numerical method 
varies inversely to x- I^^ relation to this, the number 
of operations necessary to update the MPS representa- 
tion every time a unitary transformation is applied grows 
with the size of the local basis (M -I- 1), but is attenu- 
ated by exploiting conservation of number of particles. 
Moreover, from the arguments in [ll| it follows that the 
number of operations required to update the state must 
grow as a polynomial of %. This makes Unfolding suitable 
for systems with little entanglement. However, because 
every time the state is computed only one round of uni- 
tary operations is invoked (Eq. [TT]) . Unfolding is different 
to methods where the calculation of the state for a given 
time entails an integration of short evolutions. The fact 
that in Unfolding error does not accumulate with time is 
also an advantage, as well as the fact that specific choices 
of boundary conditions or number of neighbors do not 
necessarily preclude the application of the method. The 
key point is to put the state in the form of Eq. (|15p . 
Likewise, the advantages of Unfolding over diagonaliza- 
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FIG. 3: Boson chain with N = 100, M = 100 and the same 
conditions as in Fig. [l] In this example the size of the MPS 
representation was bounded by setting x = 50. Inset. Eigen- 
values of the single-site density matrix for different times. As 
the logarithmic plot of the eigenvalue distribution becomes 
more linear, the state approaches a Gaussian state. 



tion can be appreciated by noticing that while the basis 
of H grows exponentially with N, the bases of the ma- 
trices involved in the Unfolding calculation grow linearly 
with N. In comparison to other methods that could be 
applied in the same circumstances, Unfolding would be 
suitable when a MPS description is preferred or when 
entanglement is small. 

For states with large entanglement, Unfolding can be 
used to get an estimation. This is done by setting % to a 
numerically manageable value. Fig. [3]shows some results 
obtained by fixing x ii^ ^ simulation of a relatively big 
chain. In spite of the approximation, the relaxation pro- 
file shows good agreement with theoretical assessments 
reported in [12| . 



DISCUSSION AND OUTLINE 

Although Unfolding has been presented in the context 
of a specific class of initial states, it appears the same 
strategies can in principle be applied whenever the state 
is in general given by 



N N 

f\J2ck,N'al...,Y.c,,ia]] |0), (21) 



as long as / could be expanded in Taylor series. Fur- 
thermore, coherent states like 



3Ej"ja]-"j*ai|o), 



(22) 



exhibit some compatibility with Unfolding too. In 
these cases. Folding would reduce the state to a func- 



tion of ladder operators acting on |0). The success of the 
method would then depend on the possibility of writing 
such a reduced state in MPS terms without much effort. 
The translated state can then be used as the initial con- 
dition in a simulation effectuated by, for instance, TEBD. 
As commutativity of NDMs (Eq. ([T5)) ) is assumed 
in Unfolding, non-commuting NDMs can be treated by 
adding modes that correct this anomaly. As an example, 
consider the state 



E ^^M«] E '"^.^^l |0)' E ^jm42 ^ 0. (23) 
A third mode can be introduced so that 
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Cfc.24 |0>, E^^-i42 = 0- (24) 



Up to a normalization constant, the new state can 
be folded as shown above. Once the transformation to 
MPS has been carried, the coefficients related to the ex- 
tra mode can be dropped. This approach is resembling 
of density matrix purification. On the other hand, one 
way of taking interaction effects into account is to apply 
perturbation theory, treating non-linear terms as pertur- 
bations. This would result especially effective when the 
non-linearity is local, because the MPS description is ap- 
propriate to find local mean values. Another way is to 
mimic the interaction using a mean-field approach. This 
could be realized by using the solution of the non-linear 
Gross-Pitaevskii equation as the coefficients of Eq. ^. 
One can also think of using Eq. P7)) as a variational 
ansatz, similar to the Gutzwiller ansatz. 

An alternative method has been proposed in the con- 
text of linear bosonic systems to compute physical quan- 
tum states in MPS form. The technique has been used to 
simulate an understood model and the results have been 
compared against both data produced by diagonalization 
and theoretical studies. Agreement has been satisfactory 
in every case. Aspects related to the suitability and scope 
of the technique have been analyzed and complementary 
observations have been made. 
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